TAUP-2667-2001 
February 1, 2008 



Non-linear evolution and parton distributions 
at LHC and THERA energies 

M. Lublinsky,* ^) E. Gotsman, t b) e. Levin,* ^) and U. Maor,§ 

Department of Physics 
Technion - Israel Institute of Technology 
Haifa 32000, ISRAEL 

^) HEP Department 
School of Physics and Astronomy 
Raymond and Beverly Sackler Faculty of Exact Science 
Tel Aviv University, Tel Aviv, 69978, ISRAEL 

February 1, 2008 

Abstract 

We suggest a new procedure for extrapolating the parton distributions from HERA en- 
ergies to higher energies at THERA and LHC. The procedure suggested consists of two 
steps: first, we solve the non-Unear evolution equation which includes the higher twists 
contributions, however this equation is deficient due to the low {log{\/x)) accuracy of our 
calculations. Second, we introduce a correcting function for which we write a DGLAP type 
linear evolution equation. We show that this correcting function is small in the whole kine- 
matic region and decreases at low x. The nonlinear evolution equation is solved numerically 
and first estimates for the saturation scale, as well as for the value of the gluon density at 
THERA and LHC energies are made. We show that non-linear effects lead to damping of 
the gluon density by a factor of 2 3 at x ~ 10~^. 
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1 Introduction 



The standard perturbative QCD approach to deep inelastic scattering processes is based on two 
main ideas: 

• the DGLAP evolution equations |1|] for the leading twist parton distributions; 

• the belief that if we start QCD evolution from sufficiently high ^ 2 — 4 GeV^ the 
higher twist contributions are small. 

The DGLAP evolution has two principal difficulties. First, the evolution predicts a steep growth 
of parton distributions in the region of low x which contradicts the unitarity constraints [Q. 
Hence, we can expect large corrections to the DGLAP evolution equation, at least in the region 
of low X. The second problem is of a general nature for perturbative QCD. It is well known that 
any perturbative series is a asymptotic one. Therefore, the correct estimates for the errors in the 
pQCD approach is the value of the next order term. Practically, this means that if we calculated 
the leading order gluon density distribution, say xG^'^, and the next to leading order one xG^^'^, 
then the correct estimates for the errors will be AxG = xG^^^ — xG^'^, provided that the very 
same initial conditions are used. Such errors are large, especially at low aj^. 

It has been shown over the last two decades, that the generally held view on higher twist contri- 
butions is not correct [§. Taking into account only two terms in Operator Product Expansion 
for a parton distribution (xG) we have 

xG{x,Q') = a;G^^(x,Q2) + _ ^g^t^^^q2^ ^ ^ 

with the high twist term xG^^{x, Q"^) oc xG^'^{x, Q"^) ^ in the region of low x 0. The leading 
twist function xG^^{x, Q^) sharply increases at x — > 0. Hence, we cannot conclude that the 
higher twist contributions are small in the whole kinematic region even if they are small for the 
initial value of = Q^- The scale M is introduced for dimensional reasons. 

In this paper, we suggest a remedy for these two difficulties. First, we propose solving a nonlinear 
evolution equation which takes into account the most significant higher twist contributions and 
which specifies a high energy (low x) behavior of the parton densities, which are in accordance 
with the unitarity constraints. The parton distributions which we obtain are then amended 
by adding to the solution of the nonlinear equation N , a correcting function AA^ such that 
A^ = A^ + AA^. We then write a linear DGLAP-type equation for AA^. Since AA^ in our 
approach vanishes at low x we expect this function to be small and concentrated in the region of 
moderate x. Consequently, this function should be free from all difficulties inherent in the usual 
scheme solutions of the DGLAP equation. We would like to mention here that our approach is 
very close to that of Kimber, Kwiecinski and Martin However, our treatment differs in the 
practical way that the programme is realized. 

^The widely held opinion that the NLO parton distributions describe the experimental data better is based on 
the fact that we can change the initial parton distribution drastically and achieve a better description. 
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An extrapolation of the available parton distribution to the region of lower x is a practical problem 
for LHC energies. We need to know the parton distribution both for estimates of the background 
of all interesting processes at the LHC, such as Higgs production, and for the calculation of the 
cross sections of the rare processes which the LHC is likely to measure. In this paper we will show 
that non-linear evolution will provide a considerable taming of the parton distribution in the LHC 
range of energies, by a factor of two - three in comparison with predictions based on the DGLAP 
linear evolution. We show here that the non-linear evolution considerably diminishes the value of 
the gluon density at x = 10~^ -7- 10~^ which is the energy range of the THERA electron - proton 
collider 0. Therefore, our calculations provide the first estimates for the possible collective effects 
of high parton density QCD at THERA based on the correct non-linear evolution equation. 

The paper is organized as follows. In the next section we formulate our approach and write a 
nonlinear equation for N and linear equation for AN. Section 3 is devoted to a numerical solution 
of the nonlinear equation. The following Section 4 presents some estimates of the corrections 
induced by the DGLAP kernel. The linear BFKL equation is discussed in Section 5. In the 
Section 6 we summarize our results and mention our plans for the future. 



2 A new approach to DIS 

The DGLAP equation describes the gluon radiation which leads to an increase in the number of 
partons. However, when the parton density becomes large, annihilation processes become active 
and they suppress the gluon radiation and, so, they tame the rapid increase of the parton densities 
at a new saturation scale Qs{x) [|, H, 0. The theoretical approach as well as understanding the 
physics of parton saturation is one of the most challenging problem of QCD which has stimulated 
the development of new methods 0, H, ^ ||, |^, |10|, [Tl|, |12[ , surprisingly, all approaches lead to the 
same nonlinear evolution equation. 
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The equation is written for N{rj_,x;b) = Ima^\p^i^{r±,x;b) where a'^lp^ie is the amplitude for 
elastic scattering for a dipole of size r±. The total dipole cross section is given by 
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The deep inelastic structure function F2 is related to the dipole cross section 
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where the QED wave functions ^''^ of the virtual photon are well known |ll3| , |I4|, |T5|. The meaning 
of Eq. ( p.3| ) is simple since it describes the two stages of DIS||T6[. The first stage is the decay 



of a virtual photon into a colorless dipole {qq -pair) which is described by wave function '$'^* in 
Eq. (p.3|). The second stage is the interaction of the dipole with the target (adipoic in Eq. ( p.3| )). 
This equation is the simplest manifestation of the fact that the correct degrees of freedom at high 
energies in QCD, are color dipoles |ll3 |. 

In the equation (|2l] ), the rapidity Y = —Inx and Yq = — lna;o. The ultraviolet cutoff p is needed 
to regularize the integral, but it does not appear in physical quantities. In the large Nc limit 
(number of colors) Cp = Nc/2. 

Eq. (|2.1| ) has a very simple meaning: the dipole of size xio decays in two dipoles of sizes Xi2 and 
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xo2 with the decay probability given by the wave function |\E'| = 2 °\ ■ These two dipoles then 



interact with the target. The non- linear term takes into account the Glauber corrections for such 
an interaction. 



The linear part of Eq. ( |2.1D is the BFKL equation[^], which describes the evolution of the 
multiplicity of the fixed size color dipoles with respect to the energy Y. The nonlinear term 
corresponds to a dipole splitting into two dipoles and it sums the high twist contributions. Note, 
that the linear part of Eq. ( |2.1| ) (the BFKL equation) also has higher twist contributions and vise 
versa, the main contribution of the non-linear part is to the leading twist (see Ref. for general 



and 



arguments and Ref. [0 for explicit calculations). 

Eq. ( |2.1| ) was suggested in the momentum representation by Gribov, Levin and Ryskin 
it was proved in the double log approximation of perturbative QCD by Mueller and Qiu ||], in 
Wilson Loop Operator Expansion at high energies by Balitsky [|TU|, in color dipole approach |TB| 
to high energy scattering in QCD by Kovchegov [|TT| and in the effective Lagrangian approach 
for high parton density QCD by lancu, Leonidov and McLerran (see Ref. and Refs. [^j for 
previous efforts). Therefore, it gives a reliable tool for an extrapolation of the parton distribution 
to the region of low x. 

One can see that Eq. ( p.l|) does not depend explicitly on the target^] and all such dependence 
comes from the initial condition at some initial value xq. For a target nucleus it was argued in 
Ref. ^ 



111 that the initial conditions should be taken in the Glauber form: 



A^(xoi,xo;&) = iVGM(xoi,xo; 6) 



with 



A^GM(xoi,a;;6) = 1 - exp 



O^S'^^Ol ^qDGLAP 



(x,4/x^J^(b) 



(2.4) 



(2.5) 



The equation (2^) represents the Glauber-Mueller (CM) formula which accounts for the multiple 



dipole-target interaction in the eikonal approximation [ll4| , P0| , |21| . The function S{h) is a dipole 
profile function inside the target. The value of Xq is chosen within the interval 



/ 1 N 1 

exp( ) < Xq < — — 

as 2mK 



(2.6) 



^This independence is a direct indication that the equation is correct for all targets ( hadron and nuclei ) in 
the regime of high parton density. 
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where R is the radius of the target. In this region the value of Xq is small enough to use the low x 
approximation, but the production of the gluons (color dipoles) is still suppressed as as ln(l/x) < 
1. Therefore, in this region we have the instantaneous exchange of the classical gluon fields. 
Hence, an incoming color dipole interacts separately with each nucleon in a nucleus (see Mueller 
and Kovchegov paper in Ref. p). 

For the hadron, however, we have no proof that Eq. ( |2.4| ) is correct. As far as we understand the 
only criteria in this problem (at the moment) is the correct description of the experimental data. 
We described almost all available HERA data using Eq. ( p. 41) p2t pSj , and we feel confident using 
Eq. ( |2.4|) as an initial condition for Eq. ( p.l|) . In our model we use the Gaussian S{b) = e~^^/^^ 
form for the profile function of the hadron. The parameter i? is a phenomenological input, while 
the gluon density xG^'^^^^ is a solution of the DGLAP equation. For xq we have Eq. ( |2.6|) for a 
hadron target as well, but practically we choose xq = 10~^ which satisfies Eq. ( p.6|) and for which 
we much experimental data exists to check our initial conditions. 

Solutions to the equation (|2.1| ) were studied in asymptotic limits in Ref. p4|. A first numerical 



attempt to solve an equation similar to the equation ( ^.Ij ) was reported in Ref. However, 
the solution was only obtained for nuclei, and the main emphasize was put on extremely small 
values of x. In the following section we report on our approach to find a numerical solution of 
the equation ( |2.1| ) which differs completely from the one adopted in p5|, as we use the coordinate 



representation in which the initial conditions are of a very simple form (see Eq. ( |2.4|) ). The 
second reason for using the coordinate representation is the fact that all physical observables 
can be expressed in terms of the amplitude for the dipole-target interaction in the coordinate 
representation (see below for xG). 

Unfortunately, equation ( p.l|) is an approximate one. It sums large Inx contributions only. The 
situation can be improved at small distances, as the exact x dependence of the kernel is known 
as this is the DGLAP kernel. An attempt to obtain the elastic amplitude based on elements 
of both the equation (^TT|) and DGLAP equation was presented in Ref. The authors of that 
paper first solve a generalized DGLAP-BFKL linear equation p6[, and then add to the solution 



a nonlinear perturbation of the form presented in the equation (|2.1| ). This approach actually 
incorporates the high twist contributions in the standard way, treating them as corrections to the 
leading one. 

We suggest a different approach to the problem. First, all twist contributions should be summed 
by solving equation ( ^TTI) . We denote by iV a solution of equation (^TTp . Second, we add to the 
obtained solution a correcting function AA^, which will account for the DGLAP kernel: 

N = N + AN . (2.7) 

Assuming AA^ to be small relatively compared to A^, we propose the following linear equation for 
AA^(r_L, X] h) = r\ B{2/r^, x\ b): 

dB{Q,x;b) Gpas ^ r \ u(n ^ h\ a fo q\ 

- Pg^g[z) B[Q,-]b)dz - (2.8) 



dilnQ"^) vr Jx/xo ' z'' 

IT Jx/xo Z Z Z Arc Jx/xo V zj z 
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Here Pg^g{z) stands for the usual gluon splitting function: 

"1 - 



+ 



+ zil-z) 



V12 IP ' ^ ' 



(2.9) 



Equation (|2.8| ) is a linear equation valid in the leading InQ^ approximation, with = A/r^. The 
last term in the equation represents the correction which is due to the substitution of the BFKL 
kernel 1/^ by the correct DGLAP kernel. The first term on the right hand side of the equation ( p. 81) 
is the DGLAP evolution for the correcting function AA^, while the second term is the 'nonlinear 
"interaction" of the solutions. The initial condition AA^(r^O; 3;; 6) = N{r±o,x]b) — N{r±Q, x;b) 
is a phenomenological input at some initial distance r±o = 2/Qo to be specified. In the present 
paper we report on our solution of equation (|2.1|) , while leaving the question of the determination 
of AA^ to another paper. 



3 Numerical solution of the nonlinear equation 



In this section we report on the exact numerical solution of the equation with the initial 
condition ( p.4|) . For the transverse hadron size the value = 10 (GeV~^) is taken, this cor- 
responds to the value which is obtained from "soft" high energy phenomenology |T^, ^] and is 
in agreement with HERA data on J/\l/ photo-production p8|. For xG^'^^^^{x,Q'^) we use the 



GRV'94 parameterization and the leading order solution of the DGLAP evolution equation [p9| . 
The kinematic region where the solution of is found, ranges in x from 10~^, where the initial 



conditions are set, to x = 10~ . The maximal distance is taken to be one fermi. The value of 
the ultraviolet cutoff p is 2 x 10~^ (GeV~^). The numerical solutions obtained are checked to be 
independent of this choice. In all previous studies of equation the formal dependence on the 
impact parameter b was omitted. As our first step we too will neglect the 6-dependence. Later 
we will consider the impact parameter dependence. 



3.1 Solution without 6-dependence. 

In this subsection we assume that both sides of the equation do not depend on the impact 
parameter b. In a sense, the assumption formally corresponds to the b = case, though the whole 
approach is based on the dipole picture where the impact parameter is larger than the dipole 
sizes. 

We propose to solve the equation (|2.1| ) by the method of iterations. Equation ( ^J]) can be 
rewritten 

Ar,+i(xoi,F) = Ar,+i(xoi,F;6 = 0) = Ar,(xoi,yo)e ^ V^^^ + x 

(3.10) 

/ dye ~ ^^^^ / d2x2W(2iV.(xo2,y)-iV.(xo2,y)iV,(xi2,y)) 

JYo J p ^02^12 
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Table 1: Convergence of the iterations 
for the X independent zero it- 
eration. The table shows the 
maximal relative errors in per- 
cent. 



Table 2: Convergence of the iterations 
for the GM input for zero itera- 
tion. The table shows the max- 
imal relative errors in percent. 



We stop the iterations when the relative error between two iterations |(iVj_|_i — Ni)/Ni\ -C 1 in 
the kinematic region of the interest^. Practically, the iteration is stopped when the maximal 
relative error between last two iterations is smaller than one percent, which is of the same order 
as our numerical accuracy. As a zero iteration the initial condition (|2.4| ) is taken, which is a x 
independent function. 

Convergence of the solution for the constant value as = 0.25 is shown in the Fig. (|1]). As we go to 
smaller values of x the number of iterations required for the convergence increases. For x = 10-'' 
the solution is obtained after about 25 iterations. The convergence is very slow because the zero 
iteration taken is far from the true solution. Table (|l|) presents the maximal relative errors for 
various number of performed iterations. 

In fact, the number of iterations can be dramatically decreased by guessing a better zero iteration. 
The GM formula (|2.5|) provides a reasonable choice. This formula happens to be a good approx- 
imation to the solution, thus significantly reducing the process of convergence. Taking the GM 
formula as the zero iteration, the same solution is obtained again after about a dozen iterations. 
Table (0) shows the iteration convergence. 

Equation (p?T|) as well as the BFKL equation is derived for a constant strong couphng as- Though 
some attempts were made to introduce a running coupling constant into the BFKL equation it 
has not yet been done in a self consistent way. However, the use of the running coupling constant 
in equation (|2.1|) is justified in the double logarithmic approximation. Motivated by this fact, a 



solution of equation ( |2.1| ) with running as is also obtained. 

It is interesting to compare the solutions found with two of the widely used saturation models. 
The first one is the GM formula ( p.5| ). The second one is the Golec-Biernat and Wusthoff (GW) 

^ A formal question of the iteration convergence may certainly arise. Namely, when two iterations are sufficiently 
close it does not formally mean that a limit has been reached and even if that the limit exists at all. As a counter- 
example one can consider the sum of the series 1/n, which diverges though every step decreases. We have to stress 
that we completely ignore this question. Since a solution to the equation exists and its asymptotics are known, 
we are confident about correct convergence of our procedure. 
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r(fm) r(fm) 

Figure 1: The function N is plotted versus distance (in fermi). The six curves show the convergence of 
the iterations (from below 1,5, 10, 14,18, and 25 iterations). 
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saturation model [^. In that model the effective dipole cross section cr^j^ic{x,r±) describing the 
interaction of the qq pair with a nucleon has the form: 



(^diZieix,r±) = (To NGw{x,r^) 
Ro{x) = (x/xo)^/' (GeV-i) ; = 



23.03 (mb); A = 



= [l-exp{-rl/{ARl{x)))]; (3.11) 
0.288; xo = 3.04 x 10"^ 



The cross section ctq is a dimension-full parameter which properly normalizes the dipole cross sec- 
tion. A natural comparison is between our function N and the dimensionless saturation function 
Now of the GW model. The Fig. (0) shows the comparison between our solutions, GM formula, 
and GW model. It can be seen that the correct numerical solution is obtained between the GM 
formula and GW model. At large distances GW < N < GM, while at small distances the inverse 
situation is realized GM < N < GW. We would like to remind the reader that by construction 
the GM formula is a DGLAP solution at small distances. 
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Figure 2: The comparison between the solutions N, Glauber- Mueller formula, and GW model. The 
four curves correspond to two different solutions Nag=o.25 (large dashes), Nag -running (small 
dashes), Nqm (continuous line), and Nqw of the GW model (dots). 



We would like to add a comment on the behavior of the GM formula at large distances. At fixed 
X this function has a maximum and then decreases. Such a behavior is certainly unphysical and 
is just an artifact of using the GRV parameterization which itself possesses a similar behavior. 

One of the goals of the present research is to determine the saturation scale Qs{x). Despite many 
attempts, no exact mathematical definition of Qs{x) has so far has been found. However, some 
reasonable estimates can be obtained from the function N . As can be seen from the figures (0,0) 
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N behaves in a step like manner as a function of distance: at small distances it tends to zero, while 
at large distances the saturation value one is approached. For such a step-like kind of function it 
is natural to define the saturation scale as a position where N = 1/2: 



iV(2/Q, 



X) 



1/2. 



(3.12) 



The equation (|3.12|) defines the saturation scale Qs as a function of x (Fig. (|5,a)). The saturation 
scale of the GW model defined as l/Ro{x) is also displayed in the Fig. (^a). 
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Figure 3: The saturation scale Qs is plotted as a function Igx = logiQ^x). (a) - the scale obtained 
from the equation ( 3.1^ ), the dotted line (GW) corresponds to the saturation scale of the GW 
model; (h) - the equation l{3.1f\ ) is used to determine the scale. 



We have to comment that the accepted definition of the saturation scale ( 3.12|) is not quite 
consistent with the estimates obtained from the GM formula. In the latter approach Qs{x) is 
deduced from the requirement that the gluon packing factor in a cascade should be equal to unity. 
For large values of x (10~^ - 10"'*) the saturation scale thus obtained is somewhat larger than 
that shown in the Fig. (||,a). 



3.2 6-dependence of the solution. 

In this subsection we deal with the 6-dependence of the equation ( |2.1| ). The proposed iteration 
method can in fact be applied, but in this case the computational time increases dramatically. 
Instead, we assume that solution of the equation ( |2.1| ) preserves the very same 6-dependence as 
introduced by the initial conditions ( |2.4|) : 

iV(rx,x;6) = (1 - g-'^^^'"^) ^(^)) , (3.13) 

where k is related to the 6 = solution 

K(x,rx) = - ln(l - iV(rx,x,6 = 0)). (3.14) 

Of course, the assumption requires verification. In order to verify that ( |3.13|) is a solution of 
(|2.1| ) we should check if it indeed satisfies (|2.1|) . So, we plug A^ (|3.13|) into the rhs of (p?T|) and 
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then compare it with N{lhs{ ^^ ). Fig. (D shows the comparison as a function of the impact 
parameter b. A good fit is found. The matching becomes worse in the region close to saturation. 
This happens mainly due to problematic numerical behavior of k. 
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Figure 4: The b- dependence of the solution is compared with the model dependence of the equation 
(3.15). The graphs are plotted as functions of b'^ for two values of x: x = 10"'^ and x = 10~^. 



The continuous line is the anzatz while the dashed line is N{lhs{ p^ ). 

In order to make an additional estimate of the accuracy of the anzatz ( |3.13| ), a solution of the 



equation (|2.1| ) is obtained for a single value of the impact parameter. The only simplification 
which is made in this computation is the assumption that the impact parameter b is much larger 
than the dipole sizes: 

xoi < b; Xo2 < b. (3.15) 

In this case the 6-dependence of the rhs of ( |2.1| ) is significantly simplified. This is actually the very 
same approximation that was used in the previous subsection. Figure (^ presents the comparison 
of the anzatz (|3.13 ) with numerically computed solution for b^ = 10 GeV~^. 
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Figure 5: The solution of the eg. (2.1) for 6^ = lOGeV"^ is compared with the anzatz dependence 
of the equation ( 3.1^ ). The graphs are plotted as functions of distance for two values of x: 
X = 10~^ and x = 10~^. The continuous line is the anzatz ( \3.1S ), while the dashed line is 
the exact numeric solution. 



The match is quite good at moderate x, though it becomes worse at smaller x. Note, however, 
from figures (H) and (|]) we see that the curve of the anzatz ( p.l3|) is actually squeezed by the 



two curves it was compared to. This provides an indication of the overall uncertainty of the 
approximation, which we would roughly estimate not to exceed 10%-20%. 

The above procedure cannot be claimed to be a proof of the 6-factorization. Moreover, an analytic 
derivation of the equation (|3.13| ) does not seem to be feasible. However, the solution in the form 
( 3.13 ) is clearly a satisfactory approximation of the true solution of the equation (2.1). 



We can proceed now with the evaluation of the dipole cross section 

c^dipoie(2;,'^±) = 2 / rf^6iV(r_L,x,6) . (3.16) 



Having assumed ( |3.13| ), the dipole cross section has the form 

adipoie = 27ri?2 [in(ft:) + E^{k) + ^i] . (3.17) 

In equation ( p.l7| ) 7 denotes for the Euler constant, while Ei is the exponential integral function. 
The expression (|3.17|) predicts the In k growth of the dipole cross section, which is in agreement 



with the conclusions presented in Ref. ||2^. Fig. (^) presents a comparison between adipoie from 



(|3.17| ) and the GW model (|3.11|) . From Fig. (|^) it can be seen that the ratio o"dipoie/o"d?poie grows 



with decreasing x. This is in agreement with the logarithmic growth of the dipole cross section 
0"dipoie5 while the a^J^ie is saturated at small values of x. 

We define the gluon density xG by using the Mueller formula. It relates the gluon density to the 
dipole cross section: 



2 



xG{x,Q'^) = ^ [ ^ [ ^adipoie{x',r) . (3.18) 



Fig. (^ presents the comparison between the gluon density obtained by ( p.l8| ) and the DGLAP 



input value in the GRV parameterization. At small x DGLAP predicts steep growth of the gluon 
density (power law). As was expected, the shadowing corrections suppress the growth, and the x 
growth of the function xG is logarithmic which does not contradict unitarity. 



11 




150 
xG 125 

75 
50 
25 



Q'=2(GeV') 




DGLAP 



-7 -6 -5 -4 -3 




-7 -6 -5 -4 -3 
lg(x) 



Q'=10(GeV') 




DGLAP 



-7 -6 -5 -4 -3 



1200 
1000 
xG 800 
600 
400 
200 




Q'=80 (GeV') 




-7 -6 -5 -4 -3 
lg(x) 



Figure 7: The function xG is plotted versus lg(x). The small dashes correspond to the solution with 
the running as, while the large dashes are used for the constant as = 0.25. The continuous 
line is the GRV parameterization. 
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The gluon density xG defined as in ( p.l8 ) grows with decreasing x. In a sense, this observation 



contradicts the "super-saturation" of Ref. |^ where the density is predicted to vanish in the 
large rapidity limit. However, the gluon density definition accepted in Ref. |^ differs from ours, 
and hence the functions are not compatible. 

It is worth making an additional comment about the saturation scale Qs{x). In the previous 
subsection we defined it through the equation ( |3.12 ). Here we propose an alternative definition 



similar to the one adopted in the GM formula. Namely: 

2K(x,2/g,) = 1. (3.19) 

The above definition corresponds to N{2/Qs,x,b = 0) ^ 0.4, and it predicts a somewhat larger 
saturation scale Qs{x) (Fig. (H,b)) than the equation ( p.l2|) . Comparing the four curves of both 



figures (^a) and (|^,b) an uncertainty in the Qs{x) determination can be estimated. It is about 
30%-40% through a very large range of x, though reaching 100% uncertainty at x = 10^^. 



4 DGLAP correction - consistency check 

In this section we would like to make some comments regarding the consistency of our approach. 
It was argued previously that it is necessary to add a correction term AA^ to the solution of the 
nonlinear equation (|2.1| ) which we found. In turn, the function AA^ is a solution of the evolution 
equation (|2.8|). 

Consistency of the approach requires the support of AA^ only at moderate x. The function AA^ 
should give vanishing contributions at very small x. We also expect this function to decrease 
with Q^. In order to check the above conditions some asymptotic estimates can be made without 



explicitly solving the equation (|2.8| ). 

At very small x and large distances the function iV ~ 1. For this case equation ( |2.8| ) can be 
rewritten: 

^?(*£) = ^ /; - I) ( pi2/Q, + B(Q, £)) . (4,20) 

d[mQ'^) TT Jx/xo V zj \ 4 z z ) 

The main observation is that the evolution kernel entering the equation ( |4.20| ) is actually negative. 
Hence the function AA^ is a decreasing function of Q. 



We illustrate the point in a model where the anomalous dimension has the form 



7M = ^fi-l), (4.21) 



TT \ljJ 

where the anomalous dimension is defined by the Mellin transform of the splitting function P, 



9^9- 



OLsC, 



^(^) = I dzPg^g{z)z^. (4.22) 

71 Jo 
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We represent the functions B and through their inverse Melhn transforms: 

B{Q,x) = ^. I dujx--B{Q,uj)- N{2/Q,x) = ^. I 
2m Jc 2m Jc 

The equation ( |4.2CI| ) can be now expressed in the following form: 



dujx''^N{2/Q,u;) 



dB{Q,Lj) 



(i(lng2) 

The solution of the equation ([4 .231 ) is: 

B{Q,u;) = m-'^' 



TT 



^iV(2/Q,^) + B{Q,u) 



4 JQl 



(4.23) 



(4.24) 



with as = ots^c/'^- The function C{uj) should be determined by initial conditions dX Q = Qq. 
At very small x and not too large Q the function ~ 1 and it is a slowly varying function of Q. 
Hence, it can be taken outside of the integral in ( [4.24 ) and we obtain the solution 



1 + as 

As a result for the correcting function AN we have: 



"^/iiv(2/g,.)g^+ 



AN{2/Q,x) 



as 



l + as 



—N{2/Q,x) + 



(Q2)"s • 
4C{x) 



iQ 



(4.25) 



(4.26) 



Recall that in the large limit = ^/Q^. The solution found in ( [4.26] ) is in agreement with 
assumptions made. The function AA^ decreases at small distances and is of size 0{as) compared 
to A^. The argument presented above only aimed to prove the self consistency of our approach. 



However, the complete numerical solution of the equation (|2.8| ) is still valid for phenomenological 
reasons. The initial conditions for the function AA^ should be chosen to fit the experimental data. 



5 Numerical solution of the BFKL equation 

In this section we report on the numerical solution of the BFKL equation, which is the linear 

2CFas 



part of the equation (pTT 



A^bfkl(xoi,1^) = A^BFKL(xoi,lo)exp 



TT 



P 



+ (5.27) 



Cpots 



X dy exp 

JYo 



2 Cp «s , / X, 



In 



TT 



±51 



{y-y) 



f j2 ^01 
Jp Xq2 X^ 



2Nt, 



2 " - ' BFKL 
2 



X02,y) 



We restrict our analysis to 6 = case only. The initial conditions are assumed to be proportional 
to the gluon density at Xq = 10^^: 



A^bfkl(xoi,xo) 



■xG'(xo,4/xoi). 
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(5.28) 



In order to obtain a solution of the equation ( 5.27|) the same method of iterations is apphed. 
However, for this case the convergence of the iterations is extremely slow. Dozens of iterations 
are required in order to get a final result. Fig. presents solutions of the BFKL equation for 
as = 0.25 and for running as compared with the solutions of the equation ( |2.1| ) as well as of the 
DGLAP equation. 
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Figure 8: The junction A'^bfkl plotted together with the solution of the equation ( 2A ) and the DGLAP 
equation, a - A'^bfkl for running as; b - A'^bfkl for as = 0.25; c - DGLAP solution in GRV 
parameterization; d - N for as = 0.25; f - N for running as- 

It is seen from the Fig. that the BFKL equation predicts a very steep growth of the cross 
section as x decreases. The function A^bfkl rapidly diverges from the solutions of the equation 
(|2.1| ) and of the DGLAP equation. This numerical observation is in fair agreement with previous 
theoretical considerations about the dynamics of the BFKL equation. The solution of the BFKL 
equation touches the unitarity bound. 

However, the results of BFKL dynamics is unlikely ever to be seen. The nonlinear effects ac- 
counted for in the equation (2J.) set in at very small distances thus suppressing the BFKL 
divergence. 



6 Summary 

In the present work we suggested a new approach to DIS based on first summing all twist con- 
tributions in the leading Inx approximation. Then, we wrote down a linear evolution equation 
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for the correcting function which incorporates the correct DGLAP kernel in the leading InQ^ 
approximation. 

The method of iterations was successfully applied to the solution of the nonlinear integral equation 
( |2.1| ). A key observation is the nonlinearity of the equation, and that a good choice for the zero 
iteration, insures a rapid convergence. 

The BFKL equation is solved numerically. Its solution blows up for low values of x, and this 
is the origin for the possible violation of unitarity. However, the nonlinear effects suppress the 
growth of the cross section. These effects set in at such small distances that we do not actually 
ever expect to see the BFKL dynamics. 

The nonlinear BFKL-type equation (|2.1| ) is solved numerically. The solution predicts a saturation 
of the function at large distances. The dipole cross section obtained is not saturated as a 
function of x. Due to b integration it grows logarithmically with decrease of x. Fig. (6) shows 
that the solution to the non-linear evolution equation generates a dipole cross section which is 
quite different from the predictions of the Golec-Biernat and Wusthoff model [0 . This model is 
widely used for estimating saturation effects at HERA, and our calculations demonstrate that we 
should be very careful with the conclusions based on such estimates especially at lower x. 

We demonstrated that the solution to the non-linear equation is quite different from the models 



that have been used for estimates of the saturation effects p2| , |30[] (see Fig. (2)). However, all 
these models can be used as a first iteration of the non-linear equation which leads to faster 
convergence of the numerical procedure. 

We found that the non-linear effects are very important for the extrapolations of improved knowl- 
edge on the parton distributions to higher energies such as the THERA and LHC energies. Fig. 
(7) shows that damping due to non-linear effects leads to suppression of a factor of 2 -h 3 in 
comparison with the linear DGLAP evolution. 

The solution found for the nonlinear equation was used to estimate the saturation scale Qs{x). In 
spite of considerable uncertainty in the value of the saturation scale we predict that it grows with 
decreasing x starting from about l(GeV) at x = 10~^ and reaching around lO(GeV) ai x = 10"'', 
in accordance with the theoretical expectations 0, H, 0, ||. 

We proved that the correcting function AA^ is concentrated at the moderate values of x. This 
fact demonstrates the self-consistency of our approach, as AA^ cannot be large and can be treated 
using the linear DGLAP-type of equations. The full phenomenological study of the value of AA^ 
and its x and dependence will be published elsewhere. 

We hope that our approach will be useful for the extrapolation of the HERA parton distributions 
to higher energies (lower x). Of course much further work will be required before this area of 
QCD is fully understood. 
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